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We develop a method of simulating the full quantum dynamics of multi-mode multi- 
component Bose-Einstein condensates in a trap. We use the truncated Wigner repre- 
sentation to obtain a probabilistic theory that can be sampled. This method produces 
c-number stochastic equations that can be solved using conventional computational 

CJQi methods. The technique is valid for large mode occupation numbers. We give a 
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detailed derivation of methods of functional Wigner representation appropriate for 
quantum fields. Our approach describes spatial evolution of spinor components and 
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properly accounts for nonlinear losses. 
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I. INTRODUCTION 



The Wigner representation^ - - is a convenient and effective method of simulating the dy- 
namics of bosonic quantum fields 4 , including Bose-Einstein condensates (BECs)^. It works 
best in the limit of large particle number, where third-order derivative terms in the Wigner- 
Moyal time-evolution equation can be truncated, and direct diagonalization approaches^ 
become computationally impossible. Large particle number usually implies large numbers 
of field modes with significant population, which makes two-mode variational approaches^ - — 
less accurate. This technique has been applied to a number of quantum dynamics problems 
in both quantum optical^ - — and BEC systems 1 ^ - —. A comparison of theoretical predictions 
of quantum fluctuations with experiment has generally resulted in excellent agreement, pro- 
vided the large particle number criterion is met 1 ^-. Other methods such as the positive-P 
representation 1 ^ are known to work better 17 when the truncation approximation breaks 
down; however, the truncated Wigner technique is a numerically robust and useful method 
for BEC simulations, within its domain of applicability 1 ^ 1 ^. 

The phase-space treatment of multimode problems can be simplified by working with 
functional mappings to field operators rather than to single-mode operators. This approach 
was initially introduced by Graham^i^. Later it was used in a number of works^L2i 
without formally defining the corresponding transformations or accompanying theorems. 
In order to numerically calculate the approximate evolution of the Wigner function of a 
system, one has to truncate third-order derivative terms 4 - 1 *^, and project out modes with low 
occupation numbers. This further complicates the formal description of the method. Recent 
developments in ultra-cold atomic physics mean that processes like nonlinear damping, not 
considered in detail previously, have also become important. Accordingly, much of the 
mathematical derivation of these techniques is not readily available. 

In this paper we present a formal description of the application of the resulting truncated 
Wigner representation to simulating the multi-mode dynamics of Bose-Einstein condensates 
(BECs). We successively reduce the problem in its initial form, the master equation for 
bosonic field operators, to the system of stochastic differential equations, which have sig- 
nificantly lower computational complexity. While there is a price of making the truncation 
approximation, we emphasize that this is a systematic expansion in a small parameter, 1/N, 
where N is the particle number. The purpose of this paper is to provide a more rigorous 
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proof, within the functional analysis formalism, of several identities that are used for these 
derivations. We focus especially on the problem of nonlinear damping. This is a dominant 
relaxation mechanism in BEC systems, and is often ignored or (incorrectly) approximated 
using linear loss terms. 

We derive the resulting stochastic differential equations from the functional Fokker-Planck 
equations, and show when the corresponding truncation approximations are applicable. 



II. QUANTUM FIELDS AND DYNAMICS 

In this paper we consider a C-component Bose gas in D effective dimensions. The Hamil- 
tonian for this system is expressed in terms of bosonic field creation and annihilation oper- 
ators &j(x) and tyj(x), j = 1 . . . C, which obey standard bosonic commutation relations 

[4? j ,4?'i\=8 jk 5(x'-x). (1) 

Here x G M 13 is a .D-dimensional coordinate vector, we define 4fj = ^j{x) and fy' k = ty k (x') 
for brevity (the same abbreviation will be used for all functions of coordinates), and 5(x' — x) 
is a D-dimensional Dirac delta function. 



A. Quantized Hamiltonian 

The second-quantized Hamiltonian for the system, integrated with a D— dimensional 
volume measure dx, is 

H = Jdx !^]K jk V k dx'$}$p jk (x' - x)&fi k \ , (2) 

where Ujk is the two-body scattering potential, and the single-particle Hamiltonian Kj k is 

K jk = (~^V 2 + tiujj + Vjix^j 5 jk + m jk (t). (3) 

Here m is the atomic mass, Vj is the external trapping potential for spin j, huj is the internal 
energy of spin j, and Qj k represents a time-dependent coupling that is used to rotate one 
spin projection into another. 

If we impose a momentum cutoff k c and only take into account low-energy modes, the non- 
local scattering potential Uj k (x' — x) can be replaced by the contact potential Uj k 5(x f — a?) 25 , 



giving the effective Hamiltonian 

H = fdx jfJl^A + ^J*^**} , (4) 

where and ^ are field operators in the new restricted basis of low-energy modes, which 
is described in detail in the next section. For s-wave scattering in three dimensions the 
coefficient is Ujk = 4nh 2 ajk/m, where is the scattering length. In general, the coefficient 
must be renormalized depending on the momentum cutoff 19 ^, but the change is small if 
dxi ^> djk, where dx{ is the grid step in dimension i. 



B. Master Equation 

The time-evolution of the resulting quantum density matrix with particle losses included 
can be written as a Markovian master equation^ for the system: 



dp i 
dt h 



H,p\+J2*i f dxC t [p], (5) 



i 

where I = ■ ■ ■ ,lc) is a tuple indicating the number of atoms from each component 

involved in the inelastic interaction that causes the relevant loss. Here we have introduced 
local Liouville loss terms, that describe n-body collisional losses in the Markovian approxi- 
mation: 

d [p] = 2d t pdl - 6\6 lP - pd\6 x . (6) 

The reservoir coupling operators Oi are products of local field annihilation operators: 

c 

o, =o,(*)=n ^ 

describing local fX/^Li h J -body collision losses. In other words, this assumes that Zi, . . . Zc, 
particles of internal state quantum number j = 1, . . . C all collide simultaneously within the 
volume corresponding to the inverse momentum cutoff, and are removed from the Bose gas. 

This is a minimal approach to the complicated issue of particle loss, since it assumes that 
the reservoir of "lost" particles does not interact with the original Bose gas. The accuracy 
of this approach depends on such issues as the trapping mechanism. Since, for massive 
particles, the particle number is conserved in the non-relativistic limit, "lost" particles are 
simply in a different quantum state. The assumption that these particles don't interact 



with the original Bose gas is only valid if the trap is state-selective or the collision is highly 
exothermic, such that the resulting particles are able to move rapidly away. 



C. Field operators and restricted basis 

It is always necessary in interacting quantum field theory to define a renormalization 
together with a momentum cutoff. In the case of the truncated Wigner method, this is more 
important, as the truncation approximation is only valid if M <C N, where M is the number 
of modes and N the number of particles^. We therefore wish to treat this momentum cutoff 
procedure more carefully, as it has a direct effect on the number of modes, and hence on the 
validity of the truncation. For each component j we define an orthonormal basis consisting of 
4>j,n(x), where n e Bj is a mode identifier. The orthonormality and completeness conditions 
for basis functions are, respectively, 

A 

= (8) 

n 

where the exact nature of integration area A depends on the basis set. For example, A is 
the whole space for harmonic oscillator modes, or a box for plane waves. We assume that 
the integration J dx is always performed over A. 

Standard bosonic field operators from jT]) can be decomposed as 

where single mode operators dj tn obey bosonic commutation relations, the pair j, n serving 
as a mode identifier. The cutoff mentioned in the previous section will result in operating 
with some fixed subset of each component's basis. Let ML,- C B 3 - be these subsets. Restricted 
field operators contain only modes from the subset ML,-: 

*j = Yl ^>%>"- ( 10 ) 

Formally, these field operators have the functional type \Pj G FHm^ = (K D — > Hm.), where 
Hm 3 is the Hilbert space of the restricted subset of modes. 
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Because of the restricted nature of the operator, commutation relations ([1]) no longer 
apply. The following ones should be used instead: 

= 5 jk 5 Mj (x',x), (11) 

where 5m is a restricted delta-function. The definition is given in the Appendix, in Defini- 
tion El 

III. FUNCTIONAL WIGNER REPRESENTATIONS 

In this section we introduce and obtain properties of the functional Wigner representation. 
This will use a number of definitions and results from functional analysis. The relevant 
mathematical material that is used in this section is defined and its properties derived in 
Appendix [Bl 

As a starting point, we recall that the single- mode Wigner transformation of the operator 
A is defined as 

W sm [i] = ^ J d 2 Aexp(-Aa* + A*a)Tr jAD(A,A*)}, (12) 

where the displacement operator D(\, A*) = exp(Aa^ — A*a) was first introduced by Weyl^ 8 -. 
The detailed description of the Wigner function W(a, a*) = W sm [p] can be found ir*22. In 
this section we will extend this definition to the multimode case, using a functional notation. 

In order to specify domains and ranges of discussed functions, functionals and transfor- 
mations formally, we will employ a special notation. In general, F e A — > B — > C will 
denote a function F that depends on two values of types A and B, and has a value of type 
C. Note this expression at the same time describes a function that depends on a value of 
type A, and returns a function with a type B — )■ C. The types can be nested, for example 
(A—>B)—>(C—> D) denotes a function to function mapping. 

A. Definitions of functional operators 

We introduce complex functions A (x), which play the role of the characteristic c- number 
A in the single-mode case. The important part of the definition is the functional analogue 
of the displacement operator. 
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Definition 1. Functional displacement operator Dj G F M — > M Mj 

Dj[A, A*] = exp J dx (M>) - A**,) . 

It is also convenient to define the displacement functional as: 

Definition 2. Displacement functional D G ¥ Mj — > Fm 3 — > C 

D[A, A*, **] = exp y cfe (-Atf * + A*^) . 

It can be shown that the functional displacement operator has properties similar to its 
single- mode equivalent. 

Lemma 3. 

J^Dj[A, A*] = 4-[A, A'KtfJ + ^A'*) = (*? - iA")4-[A, A*], 

*-4-[A, A*] = 4-(A, A*)(*} + ±A') = m - h^Dj^A*}. (13) 



5 A'* 2 v j 2 

Proof. Proved using Baker-Hausdorff theorem and evaluating integrals. □ 

We define a general Wigner transformation as follows: 

Definition 4. A multi-component functional Wigner transformation W G (m. d — > n^Li^M., ) — > 
rij=i — >■ C is defined as 

where Aj G Fm-, and J S 2 A = J 5 2 Ax . . . 5 2 Ac- The notation |Mj| stands for the number 
of elements in the set Mj, so |Mj| is the total number of modes in all restricted mode 
subsets. This transforms a coordinate-dependent operator A on a restricted subset of a 
Hilbert space to a functional (W[A])[*, **]. 

Next we introduce the Wigner functional, which is a special case of Wigner transforma- 
tion. 



Definition 5. The Wigner functional W G rLf=i^Mj — > K. is 

W[*,**} ee W[p] = ^js 2 A fjjD^A^^^lJ X w, 
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where xw[A, A*] is the characteristic functional 



c 



X w[A,A*] = Tr^pl[D J [A 3 ,A*]j. (14) 

The Wigner functional has two important properties analogous to the single-mode case. 
The first one is used to successively transform operator products. 

Theorem 6. For any Hilbert- Schmidt operator A, ifW[A] = (W[A])[&, then 

wfijA) = + §4) w[i], w[^U] = (** - 14) w[i], 

W[i^] = - f 4) W[i], = (** + W[i]. (15) 

Proof. The proof uses Lemma [3] given above to transform the A LJ . Z)j product inside the 
trace, together with Lemma IB .91 from the Appendix to integrate by parts (because of the 
restriction on A, the traces will be square- integrable^) , effectively moving the differentials 
to their intended places. □ 

The second property complements the first one, providing a way to obtain expectations 
of operator products given the Wigner function. Again, it requires a supplementary lemma. 

Lemma 7. For any non-negative integer r and s: 

A \ r 

(16) 



A=0 

Proof. The factor corresponding to the j-th component in the displacement operator can be 
expanded as 

exp J dx(A^] - A**,) = £ -L { (| dxA&f) ' {- J dx^' }^ . (17) 

We can swap functional derivatives with both integration and multiplication by an indepen- 
dent function, so: 



(j dxA^]\ =r^'](jdxA^]\ , (18) 



This is a familiar result for functional derivative of integrals. We note here that it is cor- 
rect given our restricted functional derivative definitions since tyj can be expanded in our 
restricted basis set, by definition ffTUl) . 
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The successive application of the differential gives us 



3 

Similarly for the other differential: 



(^7) r (/^ A A t ) r = ^(^; t r. as) 



Thus, the differentiation will eliminate all lower order terms in the expansion, and all 
higher order terms will be eliminated by setting Aj = for every j, leaving only one operator 
product with the required order. □ 

Theorem 8. For any non-negative integers rj, Sj 

n^'/i^r I ) = f s 2 * (n*?^*) sj ) w [*,*% ( 2i ) 
j=i j sym j \j=i j 

where we have used the functional integration J 5 2 ty from the Definition \B.b\ 

Proof. By definition of the Wigner functional, the right hand side in the above equation can 
be written: 

= ^W Tl ' [pflf 6 % (/ S 2 *i*7(*d SjD l A i> A ^>n) 4fAi,A*] J . (22) 
Evaluating the integral over using Lemma IB. 81 



/ = Tr jpH / ^% ((-^) rj Am ^ ] ) A 2j ' ^ 

where Am, is the delta functional from the Definition IB. 71 Integrating by parts for each 
component in turn and eliminating terms which fit Lemma IB. 101 

I = Tr I pf[ J 5 2 A,A Mj [A,] (-^7)^ (^-)" /^:A,-Ar 



(24) 



A=0 



where Am ■ is a delta functional from the Definition IB. 71 Now, recognizing the final expression 
as a part of the previous result above in Lemma [7J we immediately get the statement of the 
theorem. □ 
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IV. SPECIFIC CASES OF TRANSFORMATIONS 



In order to Wigner transform the master equation (jSj), we will need several theorems about 
transformations of specific operator products. These theorems employ the expressions for 
high-order commutators of restricted field operators, which look somewhat similar to those 
for single-mode bosonic operators, or standard field operators from^i. 

Lemma 9. Commutators for restricted field operators: 



= l5 M (x',x)(^) 



Proof. Proved by induction. 



A further generalization of these relations is 



(25) 



□ 



Lemma 10. 



<W(V, x) 



df 



~o~~m (x , X 



d^f" 



(26) 



where f(z, z*) is a function that can be expanded into the power series of z and z* . 



Proof. Let us prove the first relation; the procedure for the second one is the same. Without 
loss of generality, we assume that f(^f', iff'') can be expanded in power series of normally 
ordered operators. Using Lemma [9j 



r,s 
r,s 

Y^frsrSuix^x^y- 1 ^') 11 

r,s 

Of 



5 M {x',x) 



□ 



The simplest case is the transformation of the linear part of the Hamiltonian (J2J 
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Theorem 11. 



dx¥A> k ,A] 



dx 



5 + 

k 



5tyj~ K ' 5^* ~ J 



Proof. Proved straightforwardly using Theorem [6] and the relation 



5 



s 



- 6 jk 5 M .(x,x) T 



(27) 



(28) 
□ 



Commutators with a Laplacian inside require a somewhat special treatment, because in 
general the Laplacian acts on basis functions. For our purposes we only need one specific 
case, and, fortunately, in this case the Laplacian behaves like a constant. 



Theorem 12. 



dx[&V 2 V,A] 



dx ( -^V 2 ^ + -^-V 2 ^* ) W[A). 
1 5^ 5** 1 



Proof. Proved using Theorem [6] and Lemma IB. Ill 



(29) 
□ 



The next theorem describes the transformation of the non-linear part of the Hamilto- 
nian (HI). 



Theorem 13. 



5 



dx 



5*, 



^j^ k % + -5 Mk (x, x) (5 jk V k + 



~ (*;<M3 - \5 Uk (x, x) (5 jk % + **) 



5^ 



5 5 5 1 



5 5 5 1, 



5^-5^*5^.4 5^-5^*5^*4 



5 5 5 1 



5 5 5 1 



\I>* WL41 

5^ fc 5** 5^- 4 ~ J 5^ fc 5^* 5#* 4 J > 1 ' 
Proof. The proof is the same as in the case of Theorem [TTJ 



(30) 
□ 
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Finally, the transformation of loss terms ([6]) requires some treatment. The proof makes 
use of two auxiliary lemmas. The first one will help us move functional differentials to their 
intended places (namely, to the left). 



Lemma 14. For J 7 G Fm — >■ F and any non-negative integer a, b: 



= T G)S*<*,«>(a^»,»i (3i) 

Proof. Proved straightforwardly by induction. □ 

The second lemma gives a way to simplify sums obtained from the application of the 
previous lemma. 

Lemma 15 (Sum rearrangement). For any non-negative integer I, u: 

l min(Z— u,j) I i-max(«,«) 

E E = I>" E ( 32 ) 

Proof. Obviously, the order v = j — k of factor / can vary from (say, when j = and 
k — 0) to I (when j = I and k — 0). Therefore: 

i min(J-Mj') i 

E E ./ v h nu- ^) E E 9(v+k,k), 

j=0 k=0 v=0 k£K(l,u,v) 

where the set K is defined as 

K(l,u,v) = {k\0 < j < I AO < k < min(Z —u,j)Aj — k — v} 
= {k\k < I - v AO < k < min(/ - u, v + k)}. 

It is convenient to consider two cases separately v < u and v > u. For the former case 

Kv<u — {k\k < I — v A < k < min(/ — u, k + v ) A v < u}. 

When v < u, k < I — v < I — u < min(/ — u, k + v) is always true, and the first inequation 
is redundant: 

K v < u = {k\0 < k < min(Z — u, v + k) A v < u}. 
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Splitting into two sets to get rid of the minimum function: 



K v < u = {k\v <uAk>OA({k<l-uAl-u<v + k)V{k<v + kAl-u>v + 
= {k\v <uA0<k<l-u}. 

For the latter case: 

K v>u = {k\k < I — v AO < k < mm(l — u, k + v) A v > u] 

= {k\v >uAk>OA((k<l-vAk<l-uAl-u<k + v) 

\/(k <l-v Ak <k + v Al-u> k + v))} 
= {k\v > u AO < k < I - v}. 

Thus 

K = K v < u U K v>u 

= {k\v <uA0<k<l-u}U {k\v > u AO <k <l -v} 
= {k\0 < k < I — max(«, v)}, 

which gives us the statement of the lemma. 

Finally, the loss transformation theorem can be proved. 

Theorem 16. The Wigner transformation of loss term (0|) is 



W 

where 



dxCi[A] 



h h la / C / r \ Jc / r \ fc c 

*EE-EE n (w) (w) 

j 1= 0fci=0 j c =0k c =0 \c=l v c/ ' 



L I)i)fe =(2-(-l)^_(_i)E c ^ 

C I Z c -max(jr c ,fc c ) 

H E Q(lc,Jc,k c ,m c )5^(x,x)¥r^- mc (^ c ) L - kc 



c=l \ m c =0 



and 



Qilj.k.m) ' L) "' !/M 



2i +k+m m\j\k\(l - k - m)\{l - j - m)f 
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Proof. Proved by applying Theorem [6j expanding products using binomial theorem, using 
Lemma [T41 to move differentials to front, and applying Lemma [T5l to transform the resulting 
summations. □ 



V. WIGNER TRUNCATION AND FOKKER-PLANCK EQUATION 

Now we have all necessary tools to transform the master equation (jSJ) with the Wigner 
transformation from Definition H] to the form of a partial differential equation. 

The single-particle term ([3]) is transformed using Theorem [11] and Theorem [12] (since Kj 
is basically a sum of Laplacian operator and functions of x): 



W 



dx 



5* 3 



(36) 



where the Wigner function W = W[p\. The nonlinear term is transformed with Theorem 
(assuming U kj = U jk ): 

5 ( _ _ _ . 5m Ax, x 



[J dx^m^„,p] 



dxU. 



5^„ 



5 



5V* 



5 Mj (x,x 



■(6 jk V k + V 3 



5 5 1 



5^j 5^* 5^ k 4 



2 
5 



3 K "*■ ft 1 

5 5 1 



6Vj 5^* 5% 4 



n)w. 



(37) 
(38) 



Loss terms (J6]) are transformed with Theorem | 

Assuming that Kj k , Uj k , and Ki are real- valued, all the transformations described above 
result in a partial differential equation for W of the form 



dW 



5 



c c , 

—A* 

3 =1 3 



<M-E^-E 



C C r-o 

fe^ ****** 



£> ife + o 



W 7 ". (39) 



Terms of order higher than 2 are produced both by the nonlinear term in the Hamiltonian 
and loss terms. Such an equation could be solved without additional approximations if there 
were only orders up to 3 (which means an absence of losses)^, but in most cases all terms 
except for first- and second-order ones are truncated. With the assumption of the state 
being coherent, the condition for truncation can be shown to be 19 
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(40) 



where N is the number of atoms. This condition is equivalent to^ 

5 Mj (x,x) < |^| 2 . (41) 

The use of this Wigner truncation allows us to simplify the results of Theorem [13] and 
Theorem I 



Lemma 17. Assuming the conditions for Wigner truncation are satisfied, the result of 
Wigner transformation of the nonlinear term can be written as 



Ujk - W*^* 1 ) w - (42) 



Proof. Proved by simplifying equation (1371) under condition (jHJ). □ 

Lemma 18. Assuming the conditions for Wigner truncation are satisfied, the result of 
Wigner transformation of the loss term can be written as 

wwi - E w n d¥ n °> + E + E E < 43 ) 

71=1 " 71=1 " 71=1 p=l f 

w/iere Oi = 0,[*] = n£=i ^c c - 

Proof. The proof is basically a simplification of the result of Theorem [16] under two condi- 
tions. First, because of (HIT) we neglect all occurrences of 6m, which means setting m c = 
for every c. Second, we are dropping all terms with high order differentials, which can 
be expressed as limiting Jc + ^2 k c < 2. The only combinations of j c and k c for which 
Z(j, k) is not zero are thus {j c = 5 cn , k c = 0, n G [1, C}}, {j c = 0, k c = 5 cn , n G [1, C}} and 
{jc = o~cn,k c = S cp ,n G [l,C],p G [1,C]}. These combinations produce terms with <5/<5\I/*, 
S/5^> n and 5 2 / 5^ p 5^* n respectively. Applying these conditions one can get the statement of 
the theorem. □ 

Thus the truncated Fokker-Planck equation (FPE) is: 

\ j=l 3 j=l 3 j=l k=l 3 / 



or, in matrix form: 



dW 
~dT 



J dx (-2Re (<5* ■ A) + Tr {6**6^}) W, 
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where we define the relevant coefficients in the FPE as: 

i ( C C \ dO* 

A = -* ( E*"** + E^**i**** - E«*af7°«' ^ 

fe=l k=l J I 3 



and 



VI. STOCHASTIC DIFFERENTIAL EQUATIONS 

Direct solution of the above FPE is generally impractical, and a Monte-Carlo or sam- 
pled calculation is called for. Since the diffusion matrix is positive-definite, the truncated 
Wigner function W is a probability distribution, provided it has a positive initial distribu- 
tion. Therefore the equation can be further transformed to the equivalent set of stochastic 
differential equations in Ito form. 



A. Stochastic Evolution 

General results on such transformations are given in Appendix O, as described by The- 
orem IC.4I Application of these methods to the truncated FPE f j44l) gives immediately a 
system of SDEs in Ito 33 form: 



<m = V [Adt + BdQ] , (47) 

where 

Bfl = V^^- (48) 
Here Qi is a functional Wiener process: 

Qi = X^,»» ( 49 ) 

and Zi n are, in turn, independent complex-valued Wiener processes. Alternatively, in 
Stratonovich form the SDEs look like 



d& = V[{A-S)dt + BdQ], 
16 



(50) 



where the Stratonovich 33 term has components 



c 



dQi f a 2 o ; 



) 



6 Un (x,x). 



(51) 



n=l Z 



These equations can now be solved using conventional methods, and any required expec- 
tations symmetrically ordered operator products can be obtained from their solution using 
Theorem [HJ 



where r c and s c is some set of non- negative integers, and () paths stands for the average over 
the simulation paths. 

B. Initial states 

Initial values for the numerical integration of equations (14 7p are obtained by finding the 
Wigner transformation of the density matrix for the desired initial state, and then sampling 
the initial values according to the resulting Wigner function. As an example, consider the 
simple case with a coherent initial state. 

Theorem 19. The Wigner distribution for a multi-mode coherent state with expectation 
values 0^(0) = al , nsMis 




(52) 




(53) 



Proof. The density matrix of the state is 



p = \ a £\ n e M)(a4 0) , n e M 




(54) 



Then the characteristic function is 




(55) 
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Using the properties of the displacement operator, this can be transformed to 

Xw (a<®) = J] exp(-A>W + A n (a£>)* - i|A| 2 ). (56) 

Finally, Wigner function is 

W c {^)= J] ( f d 2 X n eM~U< - («L 0) )*) + A;K - «i 0) ) - i|A| 2 )) 

TlgM ^ ' 

|M| 

nexp(-2|a^-ai f). 

□ 

The resulting Wigner distribution is a product of independent complex-valued Gaussian 
distributions for each mode, with an expectation value equal to the expectation value of the 
mode, and variance equal to \. Therefore the initial state can be sampled as 

«n = + -^fn, (57) 

where i] n are normally distributed complex random numbers with zero mean, {rj m rj n ) = 
and (rjmVn) = $m,n or, in other words, with real components distributed independently with 
variance \. This looks like adding half a "vacuum particle" to each mode. In functional 
form this can be written as 

9 j (x,0) = *f\x,0)+ 
where \&^(ic, 0) is the "classical" ground state of the system. 

VII. CONCLUSION 

We have formally derived all the equations necessary to describe BEC interferometry ex- 
periments statistically, given a master equation written in terms of field operators. We have 
provided general equations required to use the transformation, along with its application to 
the trapped BEC case. In the latter case, the resulting SDEs can be integrated numerically 
using conventional methods, and their solutions can be used to calculate all the required 
observables. 
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Appendix A: Wirtinger differentiation 



In this paper we are using differentiation of complex functions extensively. Instead of 
the classical definition of a differential which only works for holomorphic functions, we use 
Wirtinger differentiation 34 . One can find thorough descriptions of these rules, for example, 
in 3 ^; in this section we will only outline the basics. 

Definition A.l. For a complex variable z = x + iy and a function f(z) = u(x, y) + iv(x, y) 
the Wirtinger differential is 

df(z) _ 1 fdf .df 



dz 2 \dx dy 

One can easily prove that if f(z) is holomorphic, then the above definition coincides with 
the classical differential for complex functions. Wirtinger differential obeys sum, product, 
quotient, and chain differentiation rules (the former is applied as if f(z) = f(z, z*)). 

In addition, we will need an area integration over a complex variable: 

Definition A. 2. For a complex variable z = x + iy the integral 

/OO /"OO 
/ dx dy, 
-oo J — oo 

or, in other words, this stands for a two-dimensional integral over the complex plane. 

Such integration has a property similar to a Fourier transformation in real space. 

Lemma A. 3. If X is a complex variable, then for any non-negative integers r and s: 

j d 2 aa r (a*) s exp(-Xa* + X*a) = n 2 (—^j (J^j 5(ReX)S(ImX) 

Proof. First, using known Fourier transform relations, it is easy to prove that for real x and 
v, and non-negative integer n 

dvv n exp{±2ixv) =ir{Tt/2) n 5 in \x). 

Substituting a = x + iy, expanding the a r (a*) s term using binomial theorem and employing 
the above property, one can reach the statement of the lemma. □ 

Another important property is used extensively throughout the paper. 
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Lemma A. 4. If /(A, A*) is square-integrable, then for any complex a: 

j d 2 \^ (exp(-Aa* + X*a)f(X 7 A*)) = 0, 
J d 2 \^L (exp(-Aa* + \*a)f(\, A*)) = 0. 

Proof. Square-integrability of f means lim / = and lim / = 0, so the statement of the 

ReA— - >oo ImA— >oo 

lemma can be proved by transforming to real variables and integrating. □ 



Appendix B: Functional calculus 

This section outlines the functional calculus, which is heavily used throughout the paper. 
A detailed description is given in^, and here we only provide some important definitions 
and results which are used later in the paper. In this section we will use the definitions from 
Section III C\ namely the full basis B and the restricted basis M. Given the basis, we can 
define a correspondence between functions of coordinates and their representations in mode 
space. 

Definition B.l. Let F be the space of all functions of coordinates, which consists only of 
modes from M: Fm = (M, D — > C)m (restricted functions). The composition transformation 
Cm G C M ' — > Fm creates a function from a vector of mode populations: 

Cm(ol) = (j) n a n . 

neM 

The decomposition transformation C^ 1 G F — > C' M ', correspondingly, creates a vector of 
populations out of a function: 

(C M 1 [/])n = J dxfaf, neM. 

Note that for any / G F M , C M (C^[f]) = f. 

The result of any non-linear transformation of a function / G Fm is not guaranteed to 
belong to Fm- This requires explicit projections to be used with other restricted functions. 
This also applies to the delta function of coordinates. To avoid confusion with the common 
delta function, we introduce the restricted delta function. 
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Definition B.2. The restricted delta function 5m G Fm is defined as 

S U (X',X) = ^ 0n0n- 

Note that 5^(x', x) = 5 M (x, x'). 

Any function can be projected to M using the projection transformation. 
Definition B.3. Projection transformation Vm Fm 

V M [f}(x) = (C M (C^[f]))(x) = £>„ /" dx 'Kf = [ dx'6 M (x',x)f. 
Obviously, Vm = 1- The conjugate of Vm is thus defined as 

(v M [f}(x)y= f dx'5* m {x\x)r = v M [r]{x). 



Let J~[f] G Fm — 7- F be some transformation (note that the result is not guaranteed to 
belong to the restricted basis). Because of the bijection between Fm and C^ M ', T can be 
alternatively treated as a function of a vector of complex numbers T G C' M ' — » C°°: 

T{a) = C^lFlCMiat)}]. 

Using this correspondence, we can define functional differentiation. 

Definition B.4. The functional derivative ^ G (F M -» F) ->■ (R D -> F M -> F) is defined 



as 



Note that the transformation being returned differs from the one which was taken: the 
result of the new transformation is a function of the additional variable from R D (x'). This 
variable comes from the function we are differentiating by. 

Functional derivatives behave in many ways similar to Wirtinger derivatives. A detailed 
treatment can be found in^. In particular, the following useful lemma gives us the ability 
to differentiate functionals in a similar way to common functions: 

Lemma B.5. If g(z) is a function of complex variable that can be expanded into series of 
z n (z*) m , and functional JF\f, /*] = g(f, /*), 7 G F M -> F, then 57/5 f and 57/5 f* can be 
treated as partial differentiation of the functional of two independent variables f and f*. In 
other words: 

' ^ 5 M (x', x) 9g ['J , = 5 M {x',x) d9[l ' 1 



5f df ' 5f* ■ / d f* 
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Definition B.6. Functional integration j 5 2 f £ (Fm — > F) — > C is defined as 

If the basis contains an infinite number of modes, the integral is treated as a limit |M| — > oo. 

Functional integration has the Fourier-like property analogous to Lemma IA.31 but its 
statement requires the definition of the delta functional: 

Definition B.7. For a function A £ Fm the delta functional is 

A M [A] = J] S(Re\ n )5(Im\ n ), 

neM 

where A = C^ 1 [A] . 

The delta functional has the same property as the common delta function: 
5 2 AjqA]A M [A] = j d 2 \JF(\) \\ 5{Re\ n )5{\m\ n ) 



^[A]| As0 (Bl) 



Lemma B.8 (Functional extension of Lemma [A. 3j) . For \P £ Fm and A £ Fm, and for any 

non-negative integers r and s: 



5 2 yy r (y*) s exp\^J dx (-Atf* + A*^) 
^Im, (8Y(8Y Am[a] 



5 A* J \SA / 

Proof. The proof consists of expanding functions into sums of modes and applying Lemma fA. 31 
I Ml times. □ 



Lemma B.9 (Functional extension of Lemma [A.4j) . For a square-integrable functional F 



/ 52a ^( d [ A ' A *'*'**] f [ A ' A D =0 

1 5 2 A^ (D[A, A*, % **]F[A, A*]) = 0. 



Proof. Proved by expanding integrals and differentials into modes and applying Lemma |A]4l 

□ 
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Lemma B.10. For A £ Fm and a bounded functional F 



/<((teH-*9'H F H= o 
* A M(a'(-a'H F H= o 

Proof. Proved by expanding functional integration and differentials into modes and inte- 
grating separately over each A^, using the fact that any differential of the delta function is 
zero on the infinity. □ 

In order to perform transformations of master equations, we will need a lemma that justi- 
fies the "relocation" of the Laplacian (which is a part of the kinetic term in the Hamiltonian) 
inside the functional integral. 

Lemma B.ll. If J 7 e F M ->■ F, and Vn e M, x e OA 4> n (x) = 0, then 
J dx (V^) = J d^V 2 *)^*,**] 

A A 

Proof. The proof consists of a function expansion into a mode sum and an application of 
Green's first identity. □ 

Note that the above lemma imposes an additional requirement for basis functions, but in 
practical applications it is always satisfied. For example, in a plane wave basis eigenfunctions 
are equal to zero at the border of the bounding box, and in a harmonic oscillator basis they 
are equal to zero on the infinity (which can be considered the boundary of their integration 
area). We will assume that this condition is true for any basis we work with. 



Appendix C: Functional Fokker-Planck equation 

The general approach to numerical solution of the Fokker-Planck equation is to trans- 
form it to the equivalent set of stochastic differential equations (SDEs). In the textbooks 
this transformation is defined for real variables only^, while we have functional FPE with 
complex-valued functions. 

Our starting point is the reformulation of the theorem for real-valued multivariable FPE 
from— in terms of vectors and matrices: 
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Lemma C.l (FPE-SDEs correspondence in convenient form). If z T = (z% . . . zm) is a set 

of real-valued variables, Fokker-Planck equation 



dW 



-d T z aW + -Tr {d z d T z BB T } W 



dt z 2 

is equivalent to a set of stochastic differential equations in ltd form 

dz = adt + BdZ 

and to a set of stochastic differential equations in Stratonovich form 

dz = (a - s)dt + BdZ, 
where the noise-induced (or spurious) drift vector s has elements 

k,j 

Bi being the unit vector with elements (ej)j = 5ij. Here W = W(z) is a probability dis- 
tribution, a = a(z) is a vector function, B = B(z) is a matrix function (B having size 
M x L, where L corresponds to the number of noise sources), d z = (d Zl . . . d ZM ) is a vector 
differential, and dZ is a standard L-dimensional real-valued Wiener process. 

Proof. For details see 37 , sections 3.3 and 3.4. □ 

Theorem C.2. If ot T = (ai, . . . o;m) is a set of complex-valued variables, Fokker-Planck 
equation 

^- = -dlaW - dl*a*W + Tr {d a *dlBB H } W 

is equivalent to a set of stochastic differential equations in ltd form 

da. = adt + BdZ, 

or to Stratonovich form 

da = (a - s)dt + BdZ, 

where noise-induced drift term is 

8j =Tr{B H d a *ejB} , 

and dZ = (dX + idY)/V2 is an M -dimensional complex-valued Wiener process, containing 
two real-valued L-dimensional Wiener processes dX and dY . 
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Proof. Proved straightforwardly by transforming the equation to real variables and applying 
Lemma IC.ll □ 

Theorem C.3. If a^ c \ c = 1..C are C sets of complex variables a® = (aj . . . oc^), then 
the Fokker-Planck equation 

^=-£<&.«^-£w« w > ,Hr 

c=l c=l 
m=l n=l 

is equivalent to a set of stochastic differential equations in ltd form 

da {c) = a {c) dt + B {c) dZ, c = 1..C 

or to Stratonovich form 

da® = {a® - 8 ®)dt + B^dZ, 
where noise-induced drift term is 

c 



d=l 

and dZ is an L-dimensional complex-valued Wiener process. 

Proof. Proved by joining vectors from all components into one vector and applying Theo- 
rem El □ 

Theorem C.4. For a probability distribution W[^f, &*] G — y R ; a vector of transforma- 
tions A and a matrix of transformations B the functional FPE 
dW 



j dx (-2Re (6* ■ A) + Tr {6v.S%BB H }) W 



dt 

is equivalent to the set of SDEs in ltd form 

d^ = V [Ait + BdQ] 

or in Stratonovich form 

d& = V[(A-S)dt + BdQ] 

where 

Sj = Tr {B H d**ejB} , 
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Q is a vector of functional Wiener processes: 

Ql = <t>nZl,m 

and = (Vmi, ■ ■ ■ ,Pm c ) is a vector of projection transformations. 

Proof. Proved by expanding functional derivatives and applying Theorem I C. 3 1 The diffusion 
term has to be transformed in order to conform to the theorem: 

"dx [ dx'<p' jim <pt tn j2(^YBl k) J2€^ 



J2 J dx^B^yt; J dx^BfUp (CI) 



p€B,l ' 

Grouping terms back and recognizing the definition of projection transformation, one gets 
the statement of the theorem. □ 
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